%% Check coverage of CIs & Compute empirical coverage probability

function [empirical_cov, rcv_table] = check_CI_over_nu_robust(hat_theta, hat_vartheta, rep_base_se2, rep_val_se2, nu)
% INPUT: estimates hat_theta, hat_vartheta, se^2 rep_base_se2, rep_val_se2, EB estimate nu
% OUTPUT: empirical coverage probability empirical_cov, table of robust
% critical values and related components rcv_table

    N = size(hat_theta,1);
    J = size(hat_theta,2);

    % Compute VA and VB for each group i
    inv_terms   = 1 ./ (nu + rep_base_se2);            
    w           = inv_terms ./ sum(inv_terms, 2);  
    w_sq        = w.^2;                             
    VA = (1 + sum(w_sq, 2)) .* nu;                
    VB = rep_val_se2 + sum(w_sq .* rep_base_se2, 2);
    ratio = VA./VB;

    % Compute robust cv
    robust_cv = zeros(N,1);
    for i = 1:N
        mu = ratio(i);
        robust_cv(i) = cva(mu, [], 0.2);    
    end

    % Compute the upper and lower bounds
    bar_tau = sum(w .* hat_theta, 2);
    robust_length = robust_cv .* sqrt(VB);
    original_length = 1.282 .* sqrt(VA + VB);
    lower_bound = bar_tau - robust_length;
    upper_bound = bar_tau + robust_length;

    % Check whether the reported estimates lie in the predicted CIs
    hit = (hat_vartheta >= lower_bound) & (hat_vartheta <= upper_bound);

    % Compute empirical coverage prob
    empirical_cov = sum(hit)/N;
    

    % Combine rep_base_se2 and rep_val_se2 to form a N x (J+1) matrix
    mat = [rep_base_se2, rep_val_se2];

    % Convert the variances to std errors
    mat = sqrt(mat);
    
    % Compute the ratios of nu to se
    nu_se_ratio = sqrt(nu) ./ mean(mat, 2);

    % Generate table
    rcv_table = [VA, VB, robust_cv, nu_se_ratio, robust_length, original_length];

end